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In assumed probability density function (pdf) methods of turbulent combustion, the shape of the 
scalar pdf is assumed a priori and the pdf is parametrized by its moments for which model equations 
are solved. In non-premixed flows the beta distribution has been a convenient choice to represent 
the mixture fraction in binary mixtures or a progress variable in combustion. Here the beta-pdf 
approach is extended to variable-density mixing: mixing between materials that have very large 
density differences and thus the scalar fields are active. As a consequence, new mixing phenomena 
arise due to 1) cubic non-linearities in the Navier-Stokes equation, 2) additional non-linearities in 
the molecular diffusion terms and 3) the appearance of the specific volume as a dynamical variable. 

The assumed beta-pdf approach is extended to transported pdf methods by giving the associated 
stochastic differential equation (SDE). This enables the direct computation of the scalar pdf in a 
Monte-Carlo fashion. Using the moment equations, derived from the governing SDE, we derive 
constraints on the model coefficients of the SDE that provide consistency conditions for binary 
material mixing. The beta distribution is shown to be a realizable, consistent and sufficiently general 
representation of the marginal pdf of the fluid density, an active scalar, in non-premixed variable- 
density turbulent mixing. The moment equations derived from mass conservation are compared to 
the moment equations derived from the governing SDE. This yields a series of relations between 
the non-stationary coefficients of the SDE and the mixing physics. All rigorous mathematical 
consequences of assuming a beta-pdf for the fluid mass density. 

Our treatment of this problem is general: the mixing is mathematically represented by the diver- 
gence of the velocity field which can only be specified once the problem is defined. A simple example 
of the wide range of physical problems is isobaric, isothermal, large-density binary material mixing. 
A more complex one is mixing and combustion of non-premixed reactants in which the divergence is 
related to the source terms in the energy and species conservation equations. In this paper we seek 
to describe a theoretical framework to subsequent applications. We report and document several 
rigorous mathematical results, necessary for forthcoming work that deals with the applications of 
the current results to model specification, computation and validation of binary mixing of inert 
fluids. 

Keywords: Probability density function method; Variable-density turbulence; Transition to turbulence; Ac- 
tive scalar mixing; Beta distribution 



I. INTRODUCTION 

In turbulent flows density fluctuations may arise due to 
non-uniform species concentrations, temperature or pres- 
sure. We concentrate here on the first case, where dif- 
ferences in the fluid mass density (e.g. due to mixing of 
different-density species) are very large and play a crucial 
role in the developing flow. To distinguish from classi- 
cal constant-density shear-driven turbulence, we call the 
resulting turbulent flow variable-density (VD) pressure- 
gradient-driven turbulence (PGDT) due to the active and 
important role played by the density fluctuations in the 
coupled hydrodynamics and mixing processes. 

A canonical example of PGDT is the Rayleigh- Taylor 
(RT) instability of an interface between two fluids of dif- 
ferent densities, which occurs when an external accelera- 
tion force is directed opposite to the density gradient 1-3 . 
This phenomenon is present in terrestrial examples (e.g. 
in the unstable atmospheric boundary layer), in astro- 
physics (in a collapsed core of a massive star), as well as 
in engineering (in laser-driven or electromagnetic fusion) . 
PGDT phenomena also occur in mixtures of different- 



density species accelerated by large pressure gradients, 
as in supersonic injectors, gas turbines or scramjet com- 
bustion in hypersonic vehicles. 

Many engineering combustion simulations employ ei- 
ther assumed or transported probability density function 
(pdf) methods to compute the scalar mixing fields 4 . In 
assumed pdf methods, the shape of the pdf of certain 
material fields is assumed a priori and their distribution 
is parametrized by solving for appropriate moments. In 
contrast, transported pdf methods integrate the modelled 
evolution equations of the pdf 5 . The main advantage 
of these methods, compared to moment closures, is the 
mathematically exact and closed form of the chemical 
source term, which eliminates the need for closure as- 
sumptions for such highly non-linear processes. Another 
advantage of pdf methods is that they provide a higher 
level statistical description of the turbulent fields. 

Although the beta distribution had been widely em- 
ployed before, Girimaji 6 appears to have been the first 
to rigorously discuss it model for the pdf of tur- 

bulent mixing of inert passive scalars. Due to the lack 
of experimental data and the only direct numerical sim- 
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ulation (DNS) data 7 on two-scalar mixing at the time, 
his proposal was limited to mixing of passive scalars in 
stationary homogeneous isotropic constant-density tur- 
bulence. 

The recent DNS data of Livescu & Ristorcelli 8 sug- 
gests that the beta-pdf may also be an appropriate model 
for active scalar mixing in variable- density, buoyancy- or 
equivalently, pressure-gradient-driven turbulence 9 (and 
references therein) . The analysis of Livescu and Ristor- 
celli treat the mixing of different-density fluids with an 
external acceleration force (e.g. gravity). This is also the 
driver of the Rayleigh- Taylor instability. In RT flows, 
small initial perturbations of the interface between the 
fluids rapidly become a fully turbulent flow. The simu- 
lations analysed by Livescu and Ristorcelli are of a ho- 
mogeneous RT flow in which the mean density field is 
constant; the simulations were designed to study the is- 
sues of 1) transition, 2) non-equilibrium and 3) what new 
phenomena might arise when variable-density effects are 
important. 

Mixing in Boussinesq turbulence. If the densi- 
ties of two fluids undergoing mixing are commensurate, 
the momentum equation for a buoyancy-driven flow has 
the form of a Boussinesq fluid: the density fluctuation is 
small compared to the mean density and important only 
in the body force, and the pressure gradient is the static 
head. The second moment equations for the Boussinesq 
fluid are given by Ristorcelli & Clark 10 . In the mixing of 
two pure fluids, the pdf of the fluid density proceeds from 
a double-delta configuration to a single-peaked Gaussian- 
like distribution at late time. Given symmetric initial 
conditions the pdf is symmetric at all later times; the 
skewness is zero at all times in the mixing of a (homo- 
geneous) Boussinesq fluid if the skewness is zero at the 
outset. 

Mixing in variable-density turbulence. If the 

density of the fluids arc vastly different, the phrase 
variable- density (VD) is used to distinguish from the 
Boussinesq case. In VD flows the pressure field applied 
to the very-different-density fluids gives rise to differen- 
tial fluid accelerations and a variety of inertial effects 
(added mass effects being one of them). The advection 
term in the Navicr-Stokcs equation then contains cubic 
non-linearities, resulting in several non-Boussinesq fea- 
tures. For example, the mixing process becomes highly 
asymmetric as an initially symmetric pdf rapidly devel- 
ops a sizeable skewness depending on the density differ- 
ences of the two pure fluids. These issues, first found 
by Livescu & Ristorcelli 8 are further discussed in Rcf. 11 . 
One of the consequences of these new effects is the dy- 
namic importance of the specific volume and its correla- 
tions, required to close the hydrodynamical equations 9 . 
Similar correlations appear in the equations describing 
the mixing field and its moments. One of the purposes 
of this article, and the following papers on model devel- 
opment and application 12,13 , is assessing the utility of an 
assumed beta-pdf in representing these new VD effects 
that appear as correlations with the fluctuating specific 



volume. 

Mixing in non-equilibrium turbulence. A distinc- 
tive feature of the buoyancy-driven simulations of Livescu 
& Ristorcelli 8,14 is their transient nature. The flow evo- 
lution starts from a quiescent state, transitions to fully 
developed turbulence which is then followed by a final 
period of decay. The flow is highly non-equilibrium: the 
production-to-dissipation ratio, V/e, before the transi- 
tion at the inception of the flow is on the order of a 
few hundred; at the time of fully developed turbulence 
V/e ~ 1 and in the final decay V/e ~ 0. Consequently, a 
mixing model that is based on quasi-equilibrium assump- 
tions is almost certainly in error at all stages of the flow 
evolution. Another purpose of this article is assessing the 
utility of the beta-pdf approach to representing not only 
the variable-density effects but also such non-stationary, 
non-equilibrium effects. 

A. Objectives of this article 

In this paper we explore and seek to resolve several 
issues of both modeling and theoretical nature. We 
derive the rigorous mathematical consequences of one 
simple (but reasonable) assumption in variable-density 
binary mixing: the statistical distribution of the fluid 
mass density is beta. Using this theoretical framework, 
and a coupled development on the stochastic momentum 
equation 12 , a subsequent paper 13 deals with the applica- 
tions of the current results in model computations and 
validations. Here we seek: 

1. To present a stochastic differential equation (SDE) 
for the beta distribution. 

2. To establish the possibility and the requirements of 
a beta-pdf for modeling variable-density effects. 

3. To establish the possibility and the requirements 
of a beta-pdf for modeling turbulence that has 
a highly transitional and non-equilibrium nature. 
This issue is further addressed in the sequels 12,13 . 

4. Two sets of moment equations are derived: one 
from the SDE yielding a beta distribution and one 
from exact mass conservation. These two sets of 
equations are then compared to determine what 
physical processes the parameters in the SDE are 
related to. 

5. The two sets of moment equations are then com- 
pared to investigate the following question: If the 
density pdf were a beta distribution how would one 
consistently close very different terms in the mo- 
ment equations required to describe VD mixing? 

6. To explore what additional implications of an as- 
sumed mass density pdf might yield in polytropic 
media. 

Our primary focus will be to establish and document 
rigorous mathematical results regarding the above issues. 
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B. Outline of the paper 

It seems useful to give an outline so that the mathe- 
matics does not obscure the purpose, direction and mo- 
tivation of the presentation. 

1. Sec. II lays the groundwork and background for 
the model development. Because of our interest 
in high-Schmidt-numbcr mixing, the double-delta 
"no- mix" -limit pdf is given. This limit is also of 
interest in early-time mixing of two pure fluids. 

2. Sec. Ill presents a stochastic evolution equation 
and shows that it yields a beta distribution at all 
times. Apparently, this SDE has not been given in 
the literature. 

3. From the SDE the equations for the first several 
moments are derived in Sec. IV. 

4. In Sec. V we then establish that the SDE has the 
properties that allow it to represent the fluid mass 
density and to describe material mixing. The mo- 
ment equations from the SDE are compared to the 
moment equations derived from mass conservation. 
This establishes relations between the parameters 
in the SDE and the physical mixing processes. 

5. Up to this point the equations have addressed the 
statistically homogeneous flow, in which all the pri- 
mary mixing mechanisms are seen. In Sec. VI the 
equations for inhomogencous flows are given and 
discussed. 

6. Finally, in Sec. VII conclusions are drawn and the 
utility of the development is discussed. 

7. The essential results, the model equation and its 
consistency conditions for material mixing, arc 
summarized in Appendix A. Appendix B derives 
the consequences of a polytropic equation of state. 
The symmetric beta-pdf, a special case, is investi- 
gated in Appendix C. 



II. PRELIMINARY CONSIDERATIONS 

The context of the current work is given as a prelude 
to setting up the problem. Some basic results and issues 
relating to fluid mixing physics arc summarized. 

The fluid mass density as a mixing variable. Tur- 
bulent mixing of conserved scalars is a well-explored area 
in combustion theory and atmospheric pollution model- 
ing. Both moment closures and pdf methods have been 
used to predict statistics (or the full pdf) of mass con- 
centrations and mixture fractions. Mixing models have 
been developed for these quantities, because they are of 
the main interest in their applications. In flows with 
large density variations the mass density, g, is an active 
scalar 14 , thus it is reasonable to develop a mixing model 



directly for this quantity. This is especially important 
for pressure-gradient-driven flows, as the product of the 
mass flux and the mean pressure gradient is an important 
source of turbulence and reflects the fact that Lagrangian 
particles of different-density fluids accelerate very differ- 
ently in response to pressure gradients. 

The beta-pdf in variable-density mixing of a bi- 
nary mixture. In binary mixing between fluids of large 
density differences there is a one-to-one correspondence 
between the mass concentration and the fluid density. 
Due to the fact that the scalar is active and its fluctu- 
ations are not small, the momentum equation contains 
cubic non-linearities from the product of density and the 
quadratic velocity. For this reason it is more convenient 
in variable-density turbulence, from the viewpoint of the 
moment equations, to employ the fluid density as the 
mixing variable. 

The beta distribution is a natural choice to represent 
the density in binary fluid flows since g can be expressed 
MS irl function of two dependent passive scalars, as 



Qi 



Q2 ; 



(1) 



where g± and gi denote the constant densities of the 
pure fluids and Y\ + Y% = 1. This is the same situa- 
tion as with the mixture fraction 15 in conditional moment 
closure (CMC) 16 methods and the two-scalar mixing of 
Girimaji 6 in constant-density flows. In CMC the beta 
distribution is a convenient choice for representing the 
(assumed) pdf of the mixture fraction in chemically re- 
acting turbulence. In that framework the reasons for the 
choice of the beta distribution are: 

1. The mixture fraction is a conserved, bounded scalar 
representing the state of a two-component mixture 
(Yi and 1 - Yi). 

2. The ability to represent singular delta- function 
peaks of the unmixed state is required. 

3. A continuous (even at the singularities) and always 
integrablc pdf is required. 

4. In the CMC equations the full pdf of the mixture 
fraction is required. 

5. The physical pdf in the mixing of a binary non- 
premixed flow attains a wide variety of shapes, 
which can be well- approximated by the beta dis- 
tribution. 

Similarly, in many classes of combustion problems in 
which the thermodynamic pressure can be assumed large 
and constant, and in which a progress variable, < c < 1, 
is used to describe the mixing and combustion process, 
one has 



T 



1 



(2) 
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with the mixture temperature, T. the temperature of 
species "1", Ti, and the heat-release parameter, r. Then 
one can use the beta-pdf for the density field 17 . 

Mixing in these classes of flows is done by the fluctu- 
ating velocity divergence. In the binary mixing problem, 
for example, one has for the dilatation 

<i = d' = -£>(ln q'),u, (3) 

see Rcf. 14 . For the Boussinesq fluid problem, in which 
the fluctuating density is small compared to some con- 
stant mean density, g = Qq + q' = Qq+(3c, the dilatation is 
d! = —Dc'jjj. The dilatation in all these approximations 
can be derived from the equation of state, g = g(P, T, Y) 
and continuity. It is more complicated for the combus- 
tion problem, depending on various additional approxi- 
mations that are made. For the constant thermodynamic 
pressure case one has 



d' = - ^ In Q(T, Y) + A l n Q ( T , Y) , (4) 

In this case the dilatation, which accomplishes the molec- 
ular mixing, is most generally carried as the primary 
variable. In our treatment, in order to treat the gen- 
eral variable-density problem, the fluctuating dilatation 
is carried as the primary mixing variable. 



A. The density evolution in buoyant mixing 

The time evolution of the density pdf during homoge- 
neous buoyantly-driven variable-density mixing is shown 
in Fig. 4 of Ref. 8 , and reveals several unique features: 

1. The initial double-delta pdf of the unmixed state: 
In DNS of homogeneous RT mixing the computa- 
tional domain consists of random blobs of two fluids 
with different densities at rest initially, the pdf is 
double-delta. This is the same distribution also in 
the inhomogeneous RT layer at the centreline at 
early time, representing unmixed quiescent fluids 9 . 

2. Molecular mixing with persistent extrema: With 
time the initial delta-peaks diffuse into a series of 
distributions whose variance is decreased by the in- 
creasing occurrence of density events in the interior 
of the sample space. During this time there is still 
substantial unmixed fluid and the pdf is still bi- 
modal but with two different peaks at the station- 
ary extrema of the pure fluids. Only a small frac- 
tion of the ensemble corresponds to mixed fluid. 
This behaviour of the pdf cannot be represented by 
simple deterministic relaxation models, widely used 
for passive scalar mixing, such as the interaction 
by exchange with the mean (IEM) model 18 ' 19 and 
its variants. These models relax the full ensemble, 
i.e. all scalar fluctuations, towards the mean on the 
same time-scale, without accounting for changes in 



the scalar pdf, simply by reducing the distance be- 
tween the delta peaks and leaving the shape intact. 
Relaxing only a fraction of the ensemble on possibly 
different time-scales (thereby influencing the shape 
of the evolving distribution) requires a stochastic 
model. 

3. Mixing asymmetry: Peculiar to VD mixing is the 
fact that the two different delta-function peaks do 
not decrease in magnitude at the same rate: the 
heavy fluid mixes molecularly more slowly than 
the light fluid. This is due to inertial effects and 
the ability of the heavy fluid to resist deforma- 
tion, which induces a skewness into the distribu- 
tion, most noticeable at higher initial density ra- 
tios. This is, in part, responsible for the bub- 
ble/spike asymmetry in turbulent Ray leigh- Taylor 
layers 11 . 

4. The distribution becomes asymptotically Gaussian 
at late times: At later times the initially convex- 
shaped distribution evolves into concave shapes, 
and eventually the pdf arrives at an approximate 
(clipped) Gaussian close to a fully mixed state with 
an ever-decreasing variance. 

The behaviour of the physical pdf described above can 
not be reproduced with simple Gaussian mixing models 
or by deterministic relaxation models. To capture this 
behaviour, it appears, a fundamentally different repre- 
sentation is required. 



B. The no-mix double-delta distribution 

It is useful to have some asymptotic results. The var- 
ious moment relations for the double-delta distribution 
are now derived and summarized. 

The double-delta distribution, with delta-peaks at 
gi > and g2 > 0, with the asymmetry-parameter 
< S < 1, is defined as 

@(g) = (l-S)S(g- Ql ) + S5(g-g 2 ), (5) 

where S(g) denotes the Dirac-delta function. Then the 
mean, p, is given by 

g = J Q@(e)dg = (1 - S) gi + Sg 2 , (6) 

while the central moments for n > 2 are 

-g^= J(g- g) n @(g)dg = (1 - S)( Ql - g) n + S(g 2 - g) n . 

(7) 

The mean of v — 1/g is 

v= [-®{g)dg=(l-S)v 1 +Sv 2 , (8) 
J Q 
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with v\ = 1/ 1 Q\ and i>2 — 1/^2 and the central moments 
of v, for n > 2, are 



~-v) @{g)dg = {l-S)(v 1 -v) n +S(v 2 -v) r 



3.5 



Finally, the covariance of g and v becomes 
q'v' = J(g - ~g) f - — T?J ^(£<)dg = 1 - g-u = 



(9) 



= 1 



s* + (£l + ei)s(i-s) + (i- s) 2 

\Q2 9i J 



(10) 



which shows that q'v 1 < since {Q1/Q2 + Q2/ Qi) > 2. 

Based on the above, the mean, variance, skewness and 
kurtosis of g, in the special case of g\ = and g% = 1, 
are given by 



£1 = 5', 

e* = s(i-s), 



1-25 



J 3 / 2 ^5(1- 5)' 



^2 2 5(1 - 5) 



(11) 
(12) 

(13) 
(14) 



III. THE BETA-PDF MODEL 

This section presents a SDE that yields a beta dis- 
tribution. Our intention is to develop a model to rep- 
resent the fluid density field and its moments in transi- 
tional, variable-density, pressure-gradient-driven turbu- 
lent mixing with molecular diffusion. The formulation 
will prove useful in more specialized cases, such as pas- 
sive two-scalar mixing with equal (or unequal) diffusivity, 
or asymmetric mixing with negligible molecular diffusion. 



A. The beta distribution 

The beta distribution (Fig. 1) with parameters a > 
and /3 > is 



g^jl-gf- 1 
B{a,p) ' 



(15) 



where B(a,(3) denotes the value of the Euler beta func- 
tion 

B(a,0)= [ y^Hl-vY^dy, (16) 
Jo 

which ensures that the total probability ,^{g) integrates 
to unity. 

This distribution is bounded between < g < 1 with 
a possible singular behaviour at the extremes depending 



Q = 0.1. = 0.1 
a = 0.7, P = 0.6 
a = 1.0. = 1.0 
a = 2.5./3 = 3.8 
a = 8.0. = 8.0 




FIG. 1: Beta distributions with different parameters a and 
p. If a = ft the distribution is symmetric, while ft > a gives 
positive and a > ft gives negative skewness, respectively. For 
a > 1 and ft > 1 the distribution has a single peak, whereas 
a < 1 and ft < 1 corresponds to a U-shaped beta distribution. 



on a and /3, enabling the representation of the delta- 
function peaks in the initial unmixed state. The bound- 
cdness for all values of a and (3 ensures that in a trans- 
ported pdf method no fluid particles can have densities 
lower or higher than the initial densities, which is the 
correct behaviour if compressibility effects are negligible. 

In Fig. 1 the beta distribution is shown to possess the 
correct attributes for representing the density in VD-RT 
flows. The function can evolve from an initially convex to 
a concave shape through a series of possibly skewed dis- 
tributions, incorporating mixing asymmetry due to VD 
effects, c.f. Fig. 4 in Ref. 8 . 

The first two moments of the beta distribution, defined 
by Eq. (15), are 



g^{g)dg 



Q' 



{g--gY^{g)Ag 



a/3 



(a + /3) 2 (a + /3 + 1)' 



(17) 
(18) 



As is well known, the beta distribution is fully determined 
by its mean and variance. 

Based on Eqs. (17) and (18), the two parameters, a and 
/3, can be expressed in terms of the first two moments as 



-g{l--g)-g 12 



Q 
g' 2 



1-6' 

g(l-g)-g> 2 =(l-g) 



1 



(19) 
(20) 



which shows that the knowledge of g and g' 2 determines 
all moments of the beta-pdf. Here 



6 = 1 



g(l - g)'' 



(21) 



G 
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is the commonly used mix-metric 8 . 6 — represents pure 
unmixed fluids, while 9 — 1 occurs in the fully mixed 
state. 

The above development re-iterated some of the well- 
known mathematical characteristics of the beta distribu- 
tion and qualitatively justified its use as a suitable (and 
sufficiently general) representation of the fluid density 
pdf in variable-density RT flows. 



yields the more general sample space bounds. It is worth 
keeping in mind that the "1" in the diffusion term of 
Eq. (22) thus has the unit of g. It can be seen from 
Eq. (26) that S is non-dimensional, while b and k are 
inverse time-scales. Fig. 1 with Eq. (24) shows that S is 
responsible for the skewness and the ratio 



(28) 



B. The stationary stochastic differential equation 

This section presents the stochastic equation that 
yields a beta-pdf at all times, together with its equiv- 
alent Fokkcr-Planck equation (FPE) . 

The stationary solution of the Ito SDE governing the 
Lagrangian particle property, < g* < 1, 



dg*(t) = - (S - g*)dt + V*0*(1 - Q*)dW(t), (22) 

with the Wiener process dW(t), Ref. 20 , and constant pa- 
rameters b > 0, k > and < S < 1 is the beta distri- 
bution ^(g). This can be readily seen if Eq. (15) is sub- 
stituted into the FPE, equivalent to (and derived from) 
the SDE (22), see e.g. Ref. 20 , 



dt 



d_ 

dg 



in the stationary limit, d^/dt = 0, with all constant 
parameters 



and 0=(1-S)- 



(24) 



Eq. (24) shows that the SDE (22) with the three coef- 
ficients, 6, S and k, in fact represents a two-parameter 
distribution. 

The same result can be obtained by evaluating the in- 
tegral in the stationary solution for Eq. (23) 



Kg{l - g) 6XP \ k J u =o u{l - u)' 



S- 



(25) 



where Af is a normalization constant suitable chosen so 
that & integrates to unity. 

Note that the SDE (22) is a special case of 



dQ*(t) = -[S(g 2 - gi ) + ei -Q*}dt 



(26) 



+ vMe* - Q1K02 - g*)dW{t), 



for gi < g* < g2- By setting gi = and Q2 = l, Eq. (26) 
reduces to Eq. (22). For simplicity we will use Eq. (22), 
noting that the transformation 



2(21 + QV 



controls the convexity. For S/S < 1 and (1 — S)/S < 1 the 
shape is convex, while S/S>1 and (1 — S)/5> 1 result in 
a concave shape. 



C. The non-stationary stochastic differential 
equation 

The stationary SDE (22) is now extended to non- 
stationary processes. 

Eq. (22) represents a time-homogeneous process, since 
its drift and diffusion coefficients arc not explicit func- 
tions of time 20 . In such a case, as shown above, the so- 
lution of the equation is a stationary distribution, ^{g), 
whose statistics arc time-independent and the stochastic 
process g*(t) is statistically stationary. 

In the non-stationary case the coefficients, b(t), S(t) 
and n(t), are time-dependent and 



dg*(t) 



b(t) 



(S(t) - g*)dt + s/K(t)g*(l- g*)dW(t). 



(29) 

Eq. (29) represents a beta distribution whose shape 
evolves in time. The time-dependent model coeffi- 
cients govern the evolution of the shape of the pdf 
&[g;b(t),S(t),K(t)], abbreviated as ,9{g\t). The FPE 
equivalent to Eq. (29) is formally the same as Eq. (23) 
but with time-dependent coefficients, b(t), S(t) and n(t), 
and its solution is J?(g;t). 

The above development shows that Eq. (29) is a gen- 
eralization of Eq. (22) to a transient description. Con- 
sequently, special cases of stationary distributions can 
be obtained at any instant in time by setting all three 
coefficients to constant values, in which case Eq. (29) 
will reduce to Eq. (22). The solution, i.e. the pdf, will 
converge to a stationary distribution, ^(g;t) — > J? s (g), 
whose shape is determined by the constants & s , S s and 
k s . We will use this property of Eq. (29) to investigate 
the effects of the model coefficients on the shape of the 
pdf and its statistics. It is worth emphasizing that ^ s {g) 
is a stationary distribution in the mathematical sense, to 
which the solution of the SDE (29) converges with con- 
stant coefficients. This state is not to be confused with 
the physical state of the modelled flow, which may be 
statistically non-stationary at any given time. 

It is important to appreciate that b(t), S(t) and n(t) 
must be independent of g*, but they can be a function 
of any available single-point moments of the flow. For 
instance, b[g(t), g' 2 (t), g'v'(t), g'd'(t), . . . ,t], which is ab- 
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breviated by b(t). Thus in developing suitable specifica- 
tions for b, S and k, all available one-point statistics of 
the flow can be employed. 

We showed that a temporally evolving beta-pdf can be 
represented by a Lagrangian SDE (29) with time- varying 
coefficients. This establishes the possibility of modeling 
statistically transitional, non-equilibrium flows and is a 
crucial ingredient of the joint pdf model developed in 
Refs. 12 ' 13 for the highly transitional and non-equilibrium 
RT turbulence. 



D. Stochastic equations in pdf methods 

The mathematical properties of the different types of 
SDEs in pdf methods are used to put Eq. (29) in a his- 
torical context. 

From a mathematical viewpoint, the governing equa- 
tions, employed by transported pdf methods, written as 
a general R6 diffusion process for the scalar y* , 

Ay* (t) = A(y*,t)dt + y/B(y*,t)dW(t) , (30) 
can be categorized as: 

1. Deterministic equation with B(y*,t) = 0: such as 
relaxation models for molecular mixing, including 
the rEM-family 18 > 19 , 

2. Stochastic equation with constant diffusion, 
B(y*,t) = B(t): such as the Langcvin model for 
the velocity 21 , or the log- normal model for the 
kinetic energy dissipation time-scale 22 , 

3. Stochastic equation with linear diffusion, B(y*,t) = 
B(ay* + 6, t): such as the gamma distribution 
model for the turbulence frequency 23 , and, 

4. Stochastic equation with non-linear diffusion, 
B(y*,t): such as the Fokker-Planck model for dif- 
ferential diffusion 24 ' 25 , or Eq. (29). 

The drift is linear, A(y*,t) = A(ay* + b,t), in all these 
equations. 

Deterministic equations are incapable of influencing 
the shape of the evolving pdf, a highly desired property 
of mixing models 4 . Equations with constant diffusion 
and linear drift yield a Gaussian pdf. This is suitable 
for an unbounded quantity (e.g. velocity, log- frequency) , 
but unphysical for mixing. Equations with linear drift 
and diffusion give a gamma distribution, supported on a 
semi-infinite interval. 

Eq. (29) belongs to the 4th group. The diffusion term 
in Eq. (29) is a function of the sample space, g* , and it 
effects a non-linear mapping of the Gaussian Wiener pro- 
cess. The non-linear diffusion coefficient, \J ng*(l — g*), 
allows the SDE 1) to confine the process to a bounded in- 
terval on its sample space and 2) to influence the shape 
of its solution in a wide variety of ways. This allows 
Eq. (29) to capture some fundamentally different mixing 



phenomena that simple micro- mixing models, such as the 
1EM family, cannot reproduce. 

Apparently, Eq. (29) has not appeared in the litera- 
ture in this simple yet general form. Fox's 4,24 treatment 
of differential diffusion is general enough to include the 
beta-pdf, though he works with the symmetric distribu- 
tion, a special case of our treatment, where the drift re- 
laxes to the centre of the sample space instead of S(t). 
Cai & Lin 26 also give the symmetric equation in the sta- 
tionary case. As is shown, in general Eq. (29) represents 
a non-stationary skewed beta-pdf. 

The next section examines the three SDE coefficients, 
b(t), S(t) and that together determine a(t) and f3(t), 
in more detail. 



IV. MOMENT EQUATIONS FROM THE SDE 

In order to identify possible constraints on the three 
model coefficients, b(t), S(t) and n(t), so that Eq. (29) 
can correctly represent material mixing, evolution equa- 
tions for the first few statistical moments of ,^(g;t) are 
now derived from Eq. (23). 

Balance equations for statistics are derived as follows. 
Multiplying each term in Eq. (23) with g and integrating 
over all sample space 

J g ^T dg ^ • ■ • ' (31) 

produces the governing equation for the mean density, 
~g(t). The equation for the density variance, g' 2 {t), is 
obtained by multiplying with (g — ~g) 2 and integrating 
each term. For mathematical details on evaluating the 
integrals the reader is referred to Rcf. 5 . 

For clarity, the time-dependence of the model coeffi- 
cients and the derived moments are not explicitly stated 
in the following but implied. 

A. Mean: g 

Multiplying Eq. (23) by g and integrating over the in- 
terval < g < 1 gives the governing equation for the 
time evolution of the mean, "g, of the stochastic process, 
g*{t), governed by Eq. (29) as 

^ = 1( S ~S), (32) 
showing that the stationary mean is 

Q s = S, (33) 

where the subscript s, as before, denotes the stationary 
value, d/dt = 0. Thus Eq. (33) shows that specifying 
S in any way means specifying the stationary value, ~g s , 
to which the mean of the distribution will converge since 
b > 0. 
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Examining the statistics of the stationary state, such 
as that of ~g in Eq. (33), is helpful in determining the 
effect of the coefficients independent of time: it answers 
the question "what would be the shape of the pdf for given 
b, S and k?" Since the stationary moment is an explicit 
function of only the coefficients (and only implicitly of 
time) , it clearly indicates the effects of b, S and k on the 
given moment. On the right hand side of Eq. (33), S is 
understood as S s , but for clarity, only the moment under 
consideration (such as ~g s ) is marked by s. 



B. Variance: g' 



with Cs to be determined by the model for the given 
application. Note that Cs, in general, need not be a 
constant, i.e. Cs = Cs(t). 

Eq. (38) puts a stronger monotonicity constraint on 
the variance than Eq. (37) would and assures that <5 — > 
if and only if g' 2 — > 0, since ~g is bounded. Consequently, 
the constraint Eq. (38) on 5 establishes both physical and 
mathematical consistency of the variance evolution, since 
5^0 <^ g 72 ^ 0. 

Eq. (29) has been constrained to a monotonic non- 
increasing variance, a fundamental requirement of mixing 
models in homogeneous flows. 



Multiplying Eq. (23) by (g — ~g) 2 then integrating pro- 
duces the evolution equation for the variance of the beta- 
pdf, governed by Eq. (29), 



1 dg'' 



b dt 



Sq(1 - q) - (1 + 5)g' 2 = -3(1 - g) [l - 9(1 + 5)] , 

_ (34) 
indicating that g' 2 at any time will converge to its sta- 
tionary value 



7 5 



1 + 5 



5(1 - S). 



(35) 



Alternatively, the stationary moments, ~g s , g' 2 s , can be 
obtained by direct integration of the pdf as in Eqs. (17) 
and (18) and applying the equivalence between the model 
parameters (a, (3) and (b, S, k), Eq. (24). 

DNS data of homogeneous RT mixing 14 indicates that 
as the two fluids mix, the density variance decays mono- 
tonically and approaches zero in the fully mixed state. 
Eq. (35) shows that g' 2 —> if and only if 5 — > 0, inde- 
pendent of S. (As will be shown later, in the fully mixed 
limit ^ S ^1, hence the independence requirement.) 

Since both terms on the right hand side of Eq. (34) are 
always non-negative, the monotonicity of g' 2 can only be 
ensured if 



5< 



1 



(36) 



Thus in principle, Eq. (36) could be used to constrain 8 
via 



8 = Cs 



q(1 -g)-g 



= = Cs 



1 - 



(37) 



with < Cs < 1. The above specification, however, may 
render the numerical method unstable even in the case of 
a non-fluctuating inverse time-scale, 6, since K = 6b drives 
the stochastic term in Eq. (29) and the denominator of 
Eq. (37) is close to zero in the initial unmixed state. To 
eliminate this possibility we constrain 5 via 



n'2 



Cs- 



3(1-3) 



C s (l-0), 



(38) 



C. Mean: v 

Multiplying Eq. (23) by v and integrating gives the 
evolution equation for the mean, v, of the beta-pdf as 



ldv _ f S 
bdt~{ 2 



—2 



-v' 2 ) + ( | - 5 ) v. 



1 



(39) 



Since the drift term in the SDE for g*, Eq. (29), is linear, 
the drift in the equivalent SDE governing v* = 1/g* is 
non-linear. As a consequence, in contrast to g, the mo- 
ment equations involving v comprise an infinite hierarchy 
of non-closed system of differential equations. 



D. Covariance: g'v' 

Multiplying Eq. (23) by (g — ~g)(v — v) and integrating, 
or equivalently, employing the identity, g'v' = l—~g~-v, and 
using Eqs. (32) and (39), produce the evolution equation 
for the covariance g'v' as 



E. Third moment: g' 3 

Multiplying Eq. (23) by (g — ~g) 3 and integrating result 
in the model evolution equation of the third moment of 
the beta-pdf 



1 dg' 3 
ib~dt 



5(1 -23) + l(S -q) 
5(1 -23) + l(S -3) 



l + 5)g", 



( \ + 8 



3(1 -?)(!-' 



(41) 



whose stationary value, i.e. dg' 3 /dt = 0, at any time gives 
the skewness as 



=■3/2 



2y/6(l + 6) 1-2S 
1 + 26 ' y/S(l-S)' 



(42) 
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This means that for a symmetric distribution 5 = 1/2, 
while S < 1/2 and 5 > 1/2 will result in positive and 
negative skcwness, respectively. 



TABLE I: The equations governing the statistics of the non- 
stationary beta-pdf, derived from the SDE (29) with time- 
varying coefficients, b(t), S(t) and S(t) = n(t)/b(t). 



F. Fourth moment: g' 4 

Multiplying Eq. (23) by (g — ~g) and integrating give 
the equation governing the fourth moment of the beta- 
pdf as 



/4 



f dg 
2b~dt 



36g(l ~ g)g' 2 + [S - g + 38(1 - 2g)] g' 3 

- (1 + 36)^ 
38g 2 (l-g) 2 (l- 

- (f + 3*)i* 



+ [S -g + 38(l-2g)]g' 3 
(43) 



which can be used to derive the stationary value to which 
the kurtosis will converge at any point in time 



11 



= (1 + 5) 



65(1 -25) 2 



1 + 36 (1 + 35)(1 + 2,5)5(1 -5) 



(44) 



G. The fully-mixed limit: 5 







Taking the limit 6 — > of the stationary skewness and 
kurtosis in Eqs. (42) and (44) shows that by specifying 6 
via the monotonically decreasing function of Eq. (38), at 
t — > oo the skewness will vanish, while the kurtosis will 
approach the Gaussian value of 3, independent of 5. 

The above development is consistent with Girimaji's 6 
analysis of the beta-pdf in the limit of small variance: by 
satisfying the monotonicity constraint on 6, Eq. (38), the 
SDE (29) approximates a clipped Gaussian as t —> oo. 

To summarize, in the fully-mixed limit, 6 — > 0, the 
stationary values of the first four moments, obtained from 
Eqs. (33), (35), (42) and (44) will be 



Q 



I'm 



— 3 / 2 

g u 



T fm 



fm 



fin 



s, 

0, 

= 0, 

= 3. 



(45) 
(46) 

(47) 
(48) 



These equations characterise the asymptotic shape of the 
pdf at t — > oo. 



H. The no-mix limit: S — > oo 

Mathematically, the no-mix limit corresponds to 6 — > 
oo or cquivalently 6 — > 0. Strictly speaking, Eq. (15) is 



dt 
IV 

b dt 
ldv 
b dt 



t^-=Sg(l-g)-(l + 5)e'-- 



1 dg'v' 
b dt 
1 dg 1 * 
36 dt 

1211 
2b dt 



-S 



+ v' 2 ) + ( g - 5 



(1 - g'v') + 5 



S(l-2g) + -(S-g) 



-+5)g« 



= 3dg(l - g)g' 2 + [S - g + 36(1 - 2g)] g* - (1 + 3S)g' 4 



undefined for 6 = oo, i.e. a = j3 = 0, nevertheless it is 
useful to examine 8 — > oo as an asymptotic limit. This 
results in the following stationary formulas for the first 
four moments, obtained from Eqs. (33), (35), (42) and 
(44), 





l2 nm 
s 


g' 3 


nm 


— 3 / 2 

g' z 


s 




nm 




s 



s, 

5(1 - S), 

1-2S 
y/S(l-SY 

1 

5(1 - S) 



3. 



(49) 
(50) 

(51) 
(52) 



These equations characterise the shape of the pdf inde- 
pendent of molecular diffusion. They show the effect of 
mixing asymmetry on the statistics and that even if mix- 
ing is not allowed, the distribution is capable of repre- 
senting asymmetry (and skewness) by specifying the time 
evolution of the single model parameter S(t). 



I. Summary 

The equations governing the statistics of the SDE (29) 
are summarized in Table I. It is emphasized that these 
moment equations are 1) independent of any physics and 
2) a precise result of pure mathematical nature: they 
represent the governing equations of the moments of a 
non-stationary beta distribution. 

The constraints that must be satisfied by the coeffi- 
cients, b(t), S(t) and K,(t), for mathematical and physical 
realizability of a material mixing model in homogeneous 
flows, governed by Eq. (29), are 

1. Positivity of the inverse time-scale, b(t) > 0, 

2. Boundcdncss of < S(t) < 1, and, 
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3. To ensure a non- increasing variance, Eq. (38), with 
< C s < 1. 

The no-mix (or high-Sc) limit is obtained if n(t) 3> b(t) 
at all times. 

The above development provided necessary conditions 
for mathematical and physical consistency of Eq. (29) for 
its use as a material mixing model for passive scalars. In 
the next section the model will be shown to satisfy con- 
servation of mass which establishes sufficient conditions 
for its application as a model for the stochastic fluid den- 
sity. 



V. CONSISTENCY WITH MASS 
CONSERVATION 

Up to this point, the development has been purely 
mathematical in nature; we have made no restrictions 
regarding what physical quantity g can represent, only 
that it is a beta-distributed scalar. By constructing the 
time-evolution of the model coefficients S = n/b to sat- 
isfy Eq. (38), the governing SDE (29) has been confined 
to a temporal evolution of the particle property g* whose 
variance, g' 2 , cannot increase, a fundamental requirement 
of material mixing models in statistically homogeneous 
flows. The extension to inhomogeneous flows, which in- 
troduces no new small scale terms related to molecular 
mixing, is made subsequently. 

In the following, physics is introduced by assuming 
that g* represents the fluid density, g. The continuum 
form of mass conservation is used to derive the equations 
for the moments of the density distribution. Comparisons 
between the moment equations from the SDE (a mathe- 
matical model for the evolution of a non-stationary beta- 
pdf) and the moment equations from mass conservation 
(a mathematical model expressing the physical conser- 
vation principles) is required in order to relate the SDE 
parameters, b, S and k, to physical processes of molecular 
mixing. 



inhomogeneous flows, Eulerian notation is used.) In joint 
pdf methods containing the velocity, terms originating 
from the physical process of advection appear in closed 
form; these represent mean and turbulent transport and 
production/destruction. For simplicity we choose to in- 
corporate these in the Lagrangian derivative. It is worth 
emphasizing, that this also means that the equations for 
statistics will represent the rate of change along instan- 
taneous Lagrangian paths, as defined by Eq. (54). In 
the following, we assume the existence of a velocity pdf 
model. In joint pdf methods for a set of scalars, where 
the full velocity pdf is unavailable, turbulent transport 4 
and, as will be shown, the mass flux require closure as- 
sumptions. 

If the particle position, x*, is governed by 



dx* 



v*dt, 



(55) 



the model FPE governing the Eulerian joint pdf of den- 
sity and velocity, f(g, v; x, t), can be stated as 



df 
dt 



M- 

dxi 



d_ 

do 



;(S-Q)f 



29^ 



velocity model terms. 



Mi - Q)f] 
(56) 



The mean density equation: g. The exact equation 
in homogeneous flows, derived from Eq. (53), is 



dt 



-g'd', 



(57) 



where d! — v[ i . The evolution equation for the mean 
density along a Lagrangian path, according to the joint 
pdf model, Eq. (56), is 



do 
dt 



;(S-q). 



(58) 



sdc 



The correlation of the fluctuating density and velocity 
divergence is modelled by the mean of the drift term in 
Eq. (29) as 



A. Moment equations from conservation of mass 

The starting point for the exact equations is the conser- 
vation of mass along an instantaneous Lagrangian path 
for the density, g, and specific volume, v ~ 1/g, respec- 
tively, 

—r = —gd and — = vd, (53) 
dt dt y J 

with the dilatation d = Vi j and the Lagrangian derivative 



d 







dt dt k dxk 



(54) 



As a matter of clarity and convenience, the Lagrangian 
notation will be used, as it reduces the number of terms in 
the equations. (In a subsequent section, VI. Extension to 



g'd'=-(S-g). 



(59) 



The SDE (29), representing the fluid density, is consistent 
with conservation of mass in the mean if the above holds. 

The density variance equation: g' 2 . For a sta- 
tistically homogeneous flow the exact and the pdf-model 
density variance equations are 



dp' 2 

dt 



dg' 2 
dt 



2-g-g'd' - 2g' 2 d', 
kqO- -g)-(b + k)q 12 ~. 



(60) 
(61) 



Thus the mixing terms in the exact equation are jointly 
represented by the terms in the model equation as 



2-g-g'd' - 2g' 2 d' = kq(1 -g)-(b + n)g' 2 . (62) 
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The mean specific-volume equation: v. The ex- 
act and model equations governing the mean specific vol- 
ume are, respectively, 



dv 
d± 



dv 
dt 

sdc 



v'd>, 



= \K--S 



+ V 2 ) + (--K 



(63) 
(64) 



indicating the following relation between the dissipation 
of v and the SDE parameters: 



v'd' = [ k--S 



v 2 



K \ V. 

2 



(65) 



The density-specific-volume covariance equa- 
tion: g'v'. In variable-density turbulence the quantity 
g'v' plays a primary role in the production of the mass 
flux which, in the presence of a mean pressure gradient, 
drives the turbulence 9 . Note that 1 —~g-v = g'v'. 

The exact and pdf-model evolution equations for the 
density-specific- volume covariance are 



dg'v' 

— ; — = v-g'd' — g-v'd' 
dt 



(66) 



dg'v' 



dt 













(H 




4) 


sdc 







s 



(67) 



As above, a comparison of the right hand sides of 
Eqs. (66) and (67) relates the physical mixing processes 
to the SDE parameters. 

The third density moment equation: g' 3 . The 
exact and the pdf-model governing equations for the third 
moment arc 



l d'g 73 
3~dT 

1 ch/3 

3 ~~dT 



g' 2 -g'd' -g-g' 2 d' - g' 3 d' 



K (l-2g) + -(S-g) 



(68) 

k \ g' 3 . 
(69) 



Comparing the right hand sides shows how the SDE pa- 
rameters are related to the mixing processes: 



g' 2 -g'> 



d' -g-g' 2 d' - g' 3 d' = 
K (l-2g) + ^(S-g) 



Q' 



+ k ) g 



,,/:"! 



(70) 



The fourth density moment equation: g' A . The 
exact and modelled fourth moment equations are, respec- 
tively, 



l d^ 1 
4 dt 

1 dg 7 * 

4 ~~dF 

sde 



= g' 3 -g'd' -g-g' 3 d' - g' 4 d' 



-Kg{l-g)g' 2 - -(b + '3n)g' 4 



-[b(S-g)+Ml-2g)]g> 3 . 



(71) 



(72) 



Comparing the right hand sides shows how the SDE pa- 
rameters are related to the mixing processes: 



g' 3 -g'd' -g-g' 3 d' - g' 4 d' = 

+ ^[b(S-e) + 3K(l-2Q)]^. (73) 

The above development, relating the SDE parameters 
to physical processes, are a rigorous mathematical con- 
sequence of one assumption: the fluid mass density in a 
homogeneous flow is beta-distributed. 



B. Closure relations for g' n d', n > 1 

Eqs. (59), (62), (70) and (73) indicate that assuming a 
beta-pdf for the fluid density, the moment equations are 
inter-dependent and imply a series of relations for the 
density-dilatation statistics g' n d', n > 1. These relations 
are now explicitly detailed. 

Eq. (59) indicates that 



Q 'd' = --(S-g), 
which, when substituted into Eq. (62), implies 



g' 2 d' = - 



bg(S — g) — Kg(l - g) + (b + n)g r - 



(74) 



(75) 



Then substituting both Eqs. (74) and (75) into Eq. (70) 
results in 



g' 3 d' 



b(S-g) + ^(b + n) + K(l~2g) 



n/3 



0- 



(76) 



Eqs. (74-76), provide a series of relations for the co- 
variances of the form g' n d' , establishing a number of re- 
lations between various mixing processes and the SDE 
parameters. 

Moment and SDE relations in the no-mix limit. 

It is useful to investigate the closures, g' n d', Eqs. (74- 
76), in a perturbed state about the no- mix limit, 9 = ee, 
where < eg <C 1. This results in 



~ h Q'd' e 



2 g' 2 d'\ 



-(S-g), 
S — ~g 



b g(l — g) 1 — g 



-S+(l + 8)(l-e e ), 



(77) 
(78) 



2 g' 3 d' 



b g(l - g) 



S — g 



g+{l+26)(l-2g){l-e e ) 



- [2(S -g) + g(l + S) + 28(1 - 2g)] (1 - eg). 

(79) 

If 5 — > oo and eg — > at the same rate, Eqs. (77-79) 
result in well-defined (finite) expressions for all g' n d' in 
the no-mix state. 
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C. Connecting b, S and « to physical processes 

Applying the results of setting the right hand side of 
the moment equations that come from continuity to the 
right hand sides of the moment equations that come from 
the SDE allows one to express b, S and K in terms of 
physical mixing processes. 

Solving the relations for g' n d', n < 3, Eqs. (74-76), for 
b, S and k results in 



1 1 - F 

2 

S = Q- ^e'd', 



with 



_ 2g-g'd' + 2g' 2 d' 

Ai — = . 

g' 2 



(80) 
(81) 
(82) 

(83) 



g ,s g' 2 



A 2 = 2= + = g'd' + = - = g' 2 d 



g' 2 g' 3 



A 3 = 2g l2 -g'd' - -g-g' 2 d' - g' 3 d', 



F=l-l 



9ip 



and 



q' z q' 6 



{l-2-g)g' 2 1-25 772 3/2 ' 



(84) 
(85) 

(86) 



(87) 



Eqs. (80-82) show how the first three density-velocity- 
dilatation covariances together with the first three den- 
sity moments determine the three SDE coefficients in 
Eq. (29) or, equivalently, the two parameters of the beta- 
pdf, Eq. (24). 

The unit of Ai and A 2 is that of d! , the mixing rate, 
indicating their effect on the mixing state, while A3 cor- 
responds to asymmetry. 

Dividing Eqs. (81) and (80) yields the useful result 



5 = 



1 - 



S — g 



D l2 d' 



9 6g(l-g) r g'd' 



Eq. (88) is partitioned into two parts, responsible for the 
symmetric and non-symmetric behaviour of the pdf, re- 
spectively. The first term (from Ai) is the effect of molec- 
ular mixing, c.f. Eq. (36), and the second (from A3) is the 
mixing asymmetry. Note that Eq. (36) is a model con- 
straint, while Eq. (88) is a relationship of 5 to physics if 
the underlying pdf is beta. 

The parameters, b, S and k, of the beta-SDE, rep- 
resenting the fluid mass density, have been related to 
physical mixing processes. 



The SDE parameters in the no-mix limit. It is 

insightful to investigate the relations for b, S and k, Eqs. 
(80-82) in the no-mix (or high-Sc) limit. 

Note that $ = 1 — tp is another mix-metric (similar 
to 6) with $ = in the no-mix and $ = 1 in the fully 
mixed states. Eq. (87) shows that tp is omv defined for 
a non-symmetric distribution with finite skewness (see 
also Appendix C), which can be easily seen in the no- 
mix limit: 



= VS(l -S )_ 1-25 S ^0,-,l. 
1-2S ^5(1 - S) 2' 

(89) 

Denoting small departures from the no-mix state by 
eg > and e^ > 0, i.e. 9 = eg and tp = 1 — the SDE 
coefficients in the perturbed state are given by 



6T = 



2e</,(l - £g) 



Al + £e( / 1 ' £ ^ A 3 , 



2e^(l -eg) -sg(l -e,/,) e^(l — ee) 
2(1 -£<?)(! -£4,) 



2e 4 ,(l - eg) -eg(l - £</,) 
S™ = g- ^-JdJ. 



A 2 , 



(90) 
(91) 

(92) 



Further simplifying with e^ = eg/c = e with the positive 
constant < c < 00, Eqs. (90) and (91) become 



2 - 2ce 



Ai 



c — ce 



2 — c — c£ * 1 — a 
2 - 2(1 + c)e + 2ce 2 



A, 



(2 — c)e — ce 2 



An 



(93) 
(94) 



Eq. (93) indicates that b(t), the model mixing rate, starts 
from a well-defined finite value in the unmixed state, e — > 
0, 



-Wo 



2 - c 



Ai + cA 3 , 



(95) 



while repeatedly applying L'Hopital's rule to Eq. (94) 
leads to 



and thus 



-2A 2 , 



2(2 - c)A 2 
2Ai+c(2-c)A 3 ' 



(96) 



(97) 



Assuming the same perturbation in both mix-metrics, i.e. 
c = 1, results in 



2A, 



>0,c=l 



2Ai + A 3 



(98) 



The SDE parameters have been related to physical 
mixing processes in the perturbed no-mix limit. 
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D. Summary 

This section detailed the consequences of representing 
the fluid mass density by a beta distribution in homoge- 
neous non-stationary mixing flows. In summary, 

1. The time-inhomogeneous governing equation (29) 
representing the density is consistent with conser- 
vation of mass. 

2. A series of relations for the correlations between 
the density and the velocity dilatation of the form, 
g' n d', n > 1, has been obtained. 

3. The physical meaning of the three model coeffi- 
cients, b(t), S(t) and K,(t), have been explicitly re- 
lated to the mixing physics, as reflected in the den- 
sity moments and density-dilatation covariances. 

4. If one were to design a moment closure for the first 
few moments of the pdf, one can now relate the 
mixing statistics in the various moment equations 
to each other by one consistency principle: models 
for the mixing processes in the different moment 
equation are related to each other in a unique self- 
consistent way if the underlying pdf is beta. 



VI. EXTENSION TO INHOMOGENEOUS 
FLOWS 

The SDE (29) is now extended to inhomogeneous flows. 
Beside micro-mixing, this allows the SDE to represent 
different macro- mixed states 4 . 



A. Exact inhomogeneous equations 

The governing moment equations, derived from conti- 
nuity, are given in both Lagrangian and Eulerian frame- 
works, establishing a correspondence between the pdf ap- 
proach and moment closures. 

The equations governing the density and specific vol- 
ume statistics, derived from Eq. (53), in inhomogeneous 
flows are 



dg _ - -tjj 

2 dt 



Q' 2 -d--Q-Q'd' - Q ,2 d', 



dv _-, —r-n- 

— — = v-d + v'a, 
at 

dg'v' 

— - — = v-g'd' — g-v'd', 
dt 

= -~^-d + ~^-Jd! - g-^ 2 d 7 -'g lw d 7 , (103) 



(99) 
(100) 
(101) 
(102) 



1 do' 4 

— = -g' i -d + g' 3 -g'd' - g-Q^d' - g'^d'. (104) 

A dt 



These equations are along an instantaneous Lagrangian 
path. The Lagrangian derivative, Eq. (54), can be used 
to write Eq. (99) in equivalent forms more traditional in 
moment closures as 



dg 



■ Vi-g,i +v',g',i = -g-d - g'd 1 



(105) 



dg &g(vi) gv l 

— + — = with (vi) = — , (106) 

at oxi g 



where the Favre average is denoted by ( • ) . Similarly, 
the exact variance equation (100) in equivalent Eulerian 
form is 



dg' 2 _ dg' 2 — _ dg' 2 v[ 
- +Vi— 2g-v'/-g,i 1 - 



dt 



dxi 



= -2g' 2 -d- 2g-g'd' - g' 2 d', (107) 



where the velocity fluctuation about the Favre average is 
v'l = Vi — (vi). Expanding the advection term in Eq. (101) 
leads to 



dv 



+ (v-Vi + v'v[) n = 2v-d+2v'd'. (108) 



The equation for the density-specific-volume covariance 
in the Eulerian framework is 



dQ'v' , _ dg'v' — _ — _ 



(109) 



+ g'Ujii', + v'v[g',i = v-g'd' - g-v'd'. 



It is useful to recast these equations in the nomenclature 
used in moment closure of VD turbulence. Following 
Livescu et al. 9 , Eq. (109) can be written as 

db t 1 + 6, x <9uV — , x 
— +Vib,i= — (gai), t -g-^- + 2g-v'd' , (110) 

where v'/ = —(v^) = —cii and b = —g'v'. 

Moment closures, such as Ref. 27 , solve equations for 
a, and b. In Eq. (110) the last two terms represent tur- 
bulent transport and micro-mixing, respectively. Sim- 
ilarly, aj accounts for molecular mixing and transport, 
since (gctj),j = g'd' + v' i g',i 1 and thus the mass- flux term 
in Eq. (110) is 

i + 5 



— {gai),i = g'v[-v n +v'v' i -g, i +g'v' i v' ri + v'v^g'^ 



dv'v' s 



dxi 



l - -v-g'd! -g-v'd', (111) 



representing production, turbulent transport and molec- 
ular mixing. 

In contrast, the left hand sides of the Eulerian equa- 
tions (105, 107, 108 and 109) explicitly detail the terms 
originating from the physical process of advection. In the 
pdf framework only the terms representing small scale 
mixing (g'd', v'd' and g' 2 d' on the right hand sides) re- 
quire closure assumptions. As described earlier, this is 
done consistently with the beta-pdf for the density, Eqs. 
(59, 62, 66 and 67). 
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B. Extension of the beta-SDE 

The extension of the SDE (29) to inhomogeneous flows 
is now attempted and its derived moment equations are 
compared to the ones derived from mass conservation. 

Compared to the homogeneous moment equations, 
(57), (60), (63), (68) and (71) (with the exception of gV) 
there is a single additional term in Eqs. (99-104), propor- 
tional to the mean dilatation, d. Interestingly, Eqs. (66) 
and (102) show that the exact equation governing g'v' is 
the same in both homogeneous and inhomogeneous flows. 

An extension of the SDE (29) to inhomogeneous flows 
should satisfy the following constraints: 

1. There must be a new drift term, otherwise the ho- 
mogeneous equation governing ~g is not modified. 

2. The new drift must be constant or linear in g* , 
otherwise the stationary solution is no longer beta. 

3. The drift must have the correct unit, i.e. that of g* . 

4. The drift must be tensorially correct, i.e. scalar. 

5. There may be an additional diffusion term that is 
independent of, linear, or quadratic in g* , otherwise 
the solution is no longer beta. 

A possible way to account for the large scale effects in 
the SDE (29), satisfying the above constraints, is with a 
second drift term 



Idg' 4 
4~dT 



dg*(t) 



b(x,t) 



(S(x,t) - g*) - g*d 



dt 



+ y/ K (x,t)g*(l - g*)dW(t). 



(112) 



The additional term represents the effect of mean vol- 
ume expansion and it is straightforward to show that 
Eq. (112) still results in a beta-pdf, since the drift is still 
linear in g*. Note that now the three model coefficients, 
b, S and k, are assumed to depend on the spatial coor- 
dinate, well, allowing for spatially inhomogeneous 
specifications. 

The Lagrangian moment equations, derived from the 
inhomogeneous SDE (112), become 



d| 
dt 

l d~g~v 
2 dt 

d~ 

dt 

dg'v' 



sdc 



dt 

l d^ 
3~dT 



sdc 



sdc 



— g-d + horn., (H3) 

-^2-d + hom., (114) 

v-d + horn., (H5) 

horn., (H6) 

-g^-d-~g~>2 horn., (117) 



= -g' A -d- g' 3 -g-d + hom., 



(118) 



sdc 



where only the new inhomogeneous terms are shown. 

Comparing the first four exact equations (99-102) with 
the model equations (113-116), indicates that the single 
term, — g*ddt, in the SDE creates the correct large scale 
terms, — 7j • d, —2g' 2 ■ d and v ■ d in the equations for 
7j, g' 2 and v. However, it gives rise to a spurious term, 
_ n g/(n-i) .~g . d, in the equations for g' n , n > 3, second 
terms on the right hand side of Eqs. (117, 118). There- 
fore, a consistent representation in Eq. (112) for inhomo- 
geneous flows is limited to the first two moments. 

This section extended the SDE (29) for inhomogeneous 
flows. The representation in the SDE (112) is consistent 
up to the first two density moments and for the mean 
specific volume and density-specific-volume covariance. 



VII. CONCLUSION 

The rigorous mathematical consequences of assuming 
a beta distribution for the fluid mass density in variable- 
density (VD) flows have been derived. Several issues of 
modeling and theoretical nature have been explored in or- 
der to lay the groundwork for the related and subsequent 
articles 12 ' 13 that apply these results in model computa- 
tions and validations. Our treatment have been general: 
treating the mixing process in terms of a general dilata- 
tional field. 

The main results can be summarized as follows: 

1. We have presented a stochastic differential equation 
(SDE) that yields a beta distribution and reflects 
conservation of mass in VD flows. 

2. We have derived two sets of moment equations, one 
from the SDE and one from exact mass conser- 
vation, and compared them to determine to what 
physical processes the parameters in the SDE cor- 
respond. 

3. We have drawn attention to the mean specific vol- 
ume as a primary mixing quantity, necessary to 
close the dynamical moment equations, see Ref. 9 , 
and provided a series of relations for the correla- 
tions, v'd' and g' n d' , n > 1. 

4. From a modeling point of view if one were to de- 
sign a moment closure for the first few moments of 
the pdf, one can now relate the mixing statistics in 
the various moment equations to each other using 
one consistency principle: models for the mixing 
processes in the different moment equations are re- 
lated to each other in a unique self-consistent way 
if the underlying pdf is beta. 

5. The discussion is for arbitrary Schmidt numbers. 
Some results for high Sc have been shown. 
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6. As an application, Appendix B shows that given 
the joint pdf of density (or pressure) and velocity in 
polytropic gases yields a closed diagnostic equation 
for all pressure-dilatation covariances of the form 
p' n d', n > 1. 

Our primary interest is in developing a model that can 
predict material mixing in fluids with large density dif- 
ferences, of order 10 and larger, in flows in which the tur- 
bulence is transitional and/or non-equilibrium. We have 
documented several rigorous mathematical results neces- 
sary for future works. In the subsequent papers 12 ' 13 we 
apply these ideas for VD mixing in the Rayleigh- Taylor 
instability. 
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Appendix A: Summary of results 

The solution of the stochastic differential equation, 
governing the Lagrangian particle property < g* (t) < 1 , 



dg*(t) 



-(S - q*) - g*d 



dt + y/K 8 *{l - g*)dW(t), 
(Al) 



with d = Ui^i and deterministic functions b(x,t) > 0, 
n(x,t) > and < S(x,t) < 1, is a non-stationary 
skewed beta distribution. 

For a material mixing model in turbulent flows, the 
temporal evolution of the three coefficients, b, S and k, 
needs to be specified subject to the constraint 



g(l - g) 



(A2) 



where ~g and g' 2 denote the mean and the variance, re- 
spectively, and < 9 < 1 is a mix-metric. The constraint 
(A2) ensures a non-increasing variance, as required in ho- 
mogeneous flows. 

In conjunction with a velocity pdf model, Eq. (Al) 
for the mass density is shown to be consistent with 
mass conservation up to the first two density moments 
in a variable-density, statistically inhomogencous, non- 
stationary turbulent flow. 

Special cases: 

1. Homogeneous case: The term —g*d in Eq. (Al) 
vanishes in homogeneous flows. In this case the 
model coefficients, b(t), S(t) and n(t), are only 
functions of time and Eq. (Al) is consistent with 
mass conservation for all moments of the one-point 
physical density pdf. 



2. No-mix or high- 5c limit: Restricting 5 — > oo, i.e. 
k»6, results in a model that allows mixing asym- 
metry and allows very little molecular diffusion. 



Appendix B: Example: Application to polytropic 
medium 



Throughout the paper no restriction has been made 
on the relation between the state variables: the devel- 
opment has been independent of the equation of state. 
In this section we assume that the material obeys the 
polytropic law. We show that this, together with the as- 
sumption on the density pdf, leads to some interesting 
consequences. In particular, the knowledge of the joint 
pdf of density and velocity in polytropic gases yields a 
diagnostic relationship for the pressure-dilatation covari- 
ance, p'd', which therefore requires no closure assump- 
tions, p'd' has been shown to be an important contribu- 
tor to the budget of the turbulent kinetic energy at high 
turbulent Mach numbers . 

Assuming a polytropic equation of state, 



P = Gg m , 



(Bl) 



with the polytropic coefficient, m > 0, and a constant, 
C, from mass conservation, Eq. (53), one has for the 
instantaneous pressure 



dp 
dt 



-mpd. 



(B2) 



If the fluid density is beta-distributed, an equation for 
the mean pressure, p(x, t), can be derived from the SDE 
(112) and Eq. (Bl): 



l_dp 

m dt 



sdo 



: —p-d + - [S (p-v + p'v') - p] 
+ |(m- 1) [(p-U + pV) -p] . (B3) 



Comparing this to the equation derived from mass con- 
servation, Eq. (B2), 



1 dp _ - — 

— — = -p-d-p'd', 
m dt 



(B4) 



gives a diagnostic equation for the pressure-dilatation co- 
variance as 

p'd' = — [p — S (p-v + p'v'Y\ + — (m— 1) (p — p-v — p'v') . 

(B5) 

All terms appearing in Eq. (B5) can be extracted from 
the density pdf using Eq. (Bl). For example, an equa- 
tion governing p'v' can be obtained by multiplying the 
FPE equivalent to the SDE (112) by (p -p)(v - v) and 
integrating each term. This yields 



dp'v' 



dt 



= F(b, S, k, m,p, v,p'v' ,p'v' 2 , . . . ), (B6) 



sde 
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where J 7 is a function of the SDE parameters and density 
and pressure statistics, which can all be extracted from 
the joint pdf. 

The pressure variance equations, derived from the SDE 
(112) and Eq. (B2), respectively, are 



moments, deduced from Eqs. (33), (35), (42) and (44), 



1 dp' 2 



2m dt 



sdc 



: — p' 2 -d + — (^Sp'pv — p' 
+ ^(2m - 1) (p' 2 v - p^j 



1 -ft-, 



2m dt 



p' 2 -d — p-p'd! — p' 2 d\ 



(B7) 
(B8) 



which shows a similar trend for the pressure-dilatation 
covariances, p' n d' , as for g' n d' , developed earlier: a series 
of relations can be derived between the SDE parameters, 
b, S and k, and the mixing physics but now expressed in 
terms of the pressure. The comparison of the above equa- 
tions (derived from the SDE and continuity, respectively) 
also reveal that the extension in Eq. (112) is consistent for 
the first two pressure moments in inhomogeneous flows: 
the term —g*ddt generates the correct large scale terms, 
—mp ■ d and —2mp' 2 ■ d in Eqs. (B3) and (B7). 

It is emphasized that the above development is a rig- 
orous mathematical consequence of the two assumptions: 
1) the density pdf is beta and 2) the medium is poly- 
tropic. The important point is that all the statistics g' n d' 
and p' n d' , n > 1 are known in a poly tropic medium given 
the joint pdf of density and velocity. 

To summarize: 

1. In polytropic media the joint pdf of density (or 
pressure) and velocity provides a series of consis- 
tent relations for the correlations, g' n d' and p' n d' 
for n > 1, in terms of the parameters of the pdf. 

2. Assuming a beta-pdf for the density, the equations 
governing the first two pressure moments in poly- 
tropic gases have been derived. 

3. The inhomogeneous extension of the density-SDE 
(112) consistently represents the first two moments 
of the pressure. 



Appendix C: Symmetric case: S — 1/2 

A special case of the beta-pdf can be obtained if only 
a symmetric distribution is allowed. Physically, this cor- 
responds to mixing of two equal amounts of scalars with 
equal diffusivity. 

This case is obtained by setting S = 1/2, which re- 
sults in the following stationary values for the first four 



1 

f s 



n'-'> 



— 3 / 2 

g u 



sym 



sym 



4 1 
0, 

1 - 



1/3 + 5' 



(CI) 
(C2) 

(C3) 
(C4) 



This set of equations is characteristic of the shape of the 
pdf in the symmetric case as mixing progresses. At any 
time the mean converges to 1/2, the variance decays with 
S(t), the skewness decays to zero and the kurtosis in- 
creases from 1 to 3. 

The model governing equations for the first four mo- 
ments can be obtained by setting S 1 = 1/2 in Eqs. (32), 
(34), (41) and (43), resulting in 



ldg s} 
b dt 

1 dg 12 
b ~dT 

1 



1(1 

2 V 2 



7 "(! + *) £' 2 , 



3b dt 
1 



n'3 



2b dt 



-6g' 2 - (1 + 35) £/' 



(C5) 
(C6) 

(C7) 

(C8) 



These equations show the effect of molecular diffusion in 
the special case of symmetric mixing. The time evolution 
of the single model parameter, 5(t), fully determines all 
moments of the pdf. The model coefficient b(t) acts as 
the mixing rate: the full right hand sides arc multiplied 
by multiples of b. 

Setting S = 1/2 in Eqs. (74-76) gives the closures for 
the density-velocity-dilatation covariances, g' n d', in the 
symmetric case as 



7^ ym = o, 



_ b 

~ 2 

sym b 

~ 2 



g' 2 d' 



9 ,3 d 



(l + S)g' 2 - 5 - 



^ -1(1 + S)g' 2 + (1 + 26) g' 3 



(C9) 
(CIO) 

(Cll) 



Eqs. (C9-C11) show the effect of molecular mixing on the 
statistics g' n d' , governed by the single model parameter, 
S(t), in the symmetric case. Representing molecular dif- 
fusion requires the knowledge of all the central moments, 
p 7 ™, to obtain g' n d', n>2. 
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